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Gravity  wave  characteristics  in  the  middle  atmosphere  derived 
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Abstract.  The  recently  developed  Empirical  Mode  Decomposition  (EMD)  method  is 
applied  to  analyzing  gravity  wave  characteristics  in  the  middle  atmosphere.  By  establishing 
a  close  connection  between  the  fundamental  Intrinsic  Mode  Functions  (IMFs)  derived 
from  the  EMD  method  and  WKB  solutions  of  a  dispersive-dissipative  wave  equation,  we 
show  that  the  EMD  method  can  provide  useful  insights  into  physical  processes  in  the 
middle  atmosphere  where  dispersive-dissipative  wave  phenomena  are  dominant.  A  local 
power  spectrum  function  P  is  introduced  which  provides  a  quantitative  description  of  the 
spectrum  at  any  particular  location  within  a  data  series.  The  sharp  localization  of  P  in 
space  and  wavenumber  leads  to  an  identification  of  unphysical  small  scale  oscillations  by 
falling  spheres  embedded  in  the  wind  profiles  above  60  km.  Further  analyses  of  the 
horizontal  wind  profiles  derived  from  the  Dynamics  Adapted  Network  for  the  Atmosphere 
(DYANA)  campaign  suggest  that  for  horizontal  wind  fluctuations  with  vertical 
wavenumber  m  <  3  km-1  (or  vertical  wavelength  Lz  >  2  km)  the  previously  observed 
m~ 3  Fourier  spectra  could  be  produced  by  a  linear  wave  packet  whose  characteristic 
vertical  wavenumber  decreases  with  altitude.  For  small  vertical  scale  disturbances  with 
m  >  3  km-1  ( Lz  <  2  km)  a  near  -3  slope  in  the  marginal  distribution  exists  locally  in 
the  middle  atmosphere  with  a  great  degree  of  universality,  suggesting  that  nonlinear 
energy  cascade  processes  may  dominate  the  spectral  formation  in  this  wavenumber  range. 


1.  Introduction 

In  recent  years,  internal  gravity  waves  have  been  extensively 
studied  in  the  middle  atmosphere  for  two  major  reasons.  First, 
there  are  strong  reasons  to  believe  that  gravity  waves  are  the 
primary  source  of  the  observed  mesoscale  fluctuations  in  the 
middle  atmosphere.  Second,  the  deposition  of  the  momentum 
and  energy  carried  by  the  gravity  waves  plays  a  central  role  in 
establishing  the  large-scale  circulation  and  structure  of  the 
middle  atmosphere.  Many  observations  have  suggested  a  sim¬ 
ple  form  of  m  ~p  with  p  ~  3  for  the  Fourier  vertical  wave- 
number  m  power  spectra  of  horizontal  velocity  and  relative 
temperature  fluctuations  in  the  atmosphere  at  large  m .  Many 
different  theories  have  been  developed  to  explain  the  nearly 
m  ~3  spectra  observed  at  large  m  [e.g.,  Dewan  and  Good ,  1986; 
Smith  et  al ,  1987;  Weinstock,  1990;  Hines,  1991;  Gardner ,  1994; 
Zhu ,  1994],  without  certainty  of  which  theory,  if  any,  is  valid 
[Gardner,  1996]. 

Clearly,  to  better  parameterize  the  drag  and  eddy  diffusion 
induced  by  gravity  waves  in  the  middle  atmosphere,  it  is  un¬ 


applied  Physics  Laboratory,  Johns  Hopkins  University,  Laurel, 
Maryland. 

department  of  Earth  and  Planetary  Sciences,  Johns  Hopkins  Uni¬ 
versity,  Baltimore,  Maryland. 

Computational  Physics,  Inc.,  Fairfax,  Virginia. 

4E.  O.  Hulburt  Center  for  Space  Research,  Naval  Research  Labo¬ 
ratory,  Washington,  D.  C. 

department  of  Physics,  University  of  Wuppertal,  Wuppertal,  Ger¬ 
many. 

department  of  Geophysics,  Faculty  of  Science,  Kyoto  University, 
Kyoto,  Japan. 

Copyright  1997  by  the  American  Geophysical  Union. 

Paper  number  97JD01139. 

0148-0227/97/97 JD-01 139$09.00 


portant  to  identify  the  major  physical  mechanisms  driving  the 
observed  m-3  power  spectra.  Some  of  the  difficulties  involved 
in  identifying  the  dominant  physical  mechanism  could  be  as¬ 
sociated  with  the  strong  inhomogeneity  of  the  atmosphere. 
Vertically  propagating  gravity  waves  subject  to  atmospheric 
stratification  and  wind  shear  have  nonstationaiy  wave  ampli¬ 
tudes  [e.g.,  Eckermann,  1990;  Hines,  1991].  On  the  other  hand, 
traditional  Fourier  spectral  analysis  decomposes  a  series  of 
data  into  a  linear  combination  of  infinite  harmonic  functions 
which  have  constant  amplitudes  over  the  whole  data  series. 
Thus  Fourier  methods  are  not  ideal  for  studying  nonstationary 
gravity  wave  process.  There  are  many  useful  tools  in  data 
analysis  that  provide  both  spatial  and  wavenumber  information 
to  study  nonstationary  processes  [e.g.,  Daubechies,  1992; 
Cohen,  1995].  For  example,  wavelet  analysis  has  been  applied 
to  studying  spatial  and  temporal  variations  in  gravity  wave 
activity  measured  in  the  lower  stratosphere  [Sato  and  Yamada, 
1994;  Shimomai  et  al.,  1996;  Bacmeister  et  al.,  1996,  1997]. 

Recently,  a  new  technique  called  Empirical  Mode  Decom¬ 
position  (EMD)  has  been  developed  to  study  the  nonlinear 
and  nonstationary  properties  of  time  series  [Huang  et  al.,  1997]. 
The  EMD  method  has  been  shown  to  be  superior  in  some 
circumstances  to  existing  data  analysis  tools  such  as  the  wavelet 
analysis  [Huang  et  al.,  1997].  In  this  paper,  the  EMD  method  is 
used  to  assess  the  predictions  of  current  linear  and  nonlinear 
gravity  wave  theories  of  the  observed  m~3  power  spectra  in  the 
middle  atmosphere. 

In  the  next  section  we  first  review  the  EMD  method  and 
show  that  the  fundamental  Intrinsic  Mode  Functions  (IMFs) 
derived  from  the  EMD  method  are  closely  related  to  the  WKB 
solutions  of  an  equation  describing  dispersive -dissipative  wave 
phenomena.  In  section  3,  a  two-dimensional  distribution  func¬ 
tion  called  local  power  spectrum  denoted  by  P(zk,  m€)  is 
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introduced  to  provide  a  quantitative  description  of  space- 
wavenumber  ( z-m )  localization,  where  z  and  m  are  the  alti¬ 
tude  and  vertical  wavenumber,  respectively.  Similarities  and 
differences  between  the  marginal  distribution  of  P(zk,  m€) 
and  the  Fourier  spectrum  of  a  data  series  are  clarified.  Section 
4  discusses  the  local  power  spectrum  calculated  from  the  ob¬ 
servations.  Section  5  identifies  one  IMF  as  a  characteristic 
gravity  wave  component  in  the  middle  atmosphere.  Finally,  in 
section  6  we  summarize  our  EMD-derived  results  for  gravity 
waves. 


2.  Local  Wavenumber  and  Empirical 
Mode  Decomposition  Method 

Given  a  digital  signal  or  a  data  series  of  a  quantity  u  ex¬ 
pressed  as  a  real  function  of  a  time  or  space  variable  £,  w(£), 
one  can  construct  a  complex  signal,  £/(£),  by  adding  an  imag¬ 
inary  function  /«(£)  to  the  original  signal 


U(&  =  «(f)  +  m&  =  UoU)  exp  [f*(0]. 

(1) 

where 

Udi&  =  M£)  +  «2(£)]1/2, 

(2) 

and 

Hi) 

<Ki)  =  arctan 

(3) 

are  the  amplitude  and  the  phase  of  the  signal,  respectively 
[Cohen,  1995].  When  the  independent  variable  £  represents 
space  (time),  a  local  wavenumber  (instantaneous  frequency) 
can  be  derived  from  (3) 

if/„  «(|)<}'(g)-«(f)M'(g) 

* <{)  ‘  vim 

(4) 

is  to  put  the  slowly  varying  content  of  w(£)  into  the  amplitude 
£/0(f)  and  to  put  the  fast  varying  content  into  exp  [/<£(£)]• 
Mathematically,  as  long  as  the  Hilbert  transform  (6)  exists,  one 
can  always  find  an  amplitude  U0(€)  and  a  wavenumber  k(£) 
from  (2)  and  (4)  for  any  real  function  w(£).  However,  the 
derived  f/0(f)  and  k(£)  based  on  such  a  direct  procedure  are 
not  always  physically  meaningful.  The  difficulties  are  summa¬ 
rized  by  Cohen  [1995]  into  five  paradoxes  about  the  properties 
of  the  local  wavenumber.  It  turns  out  that  the  major  difficulties 
in  finding  physically  meaningful  local  wavenumber  and  ampli¬ 
tude  for  a  real  signal  depend  on  how  well  separated  the  spectra 
of  U0( £)  and  exp  [/<#>(£)]  are-  The  recently  developed  Empir¬ 
ical  Mode  Decomposition  (EMD)  method  helps  to  resolve  this 
key  issue  by  decomposing  a  real  signal  w(f)  into  a  set  of 
different  components  so  that  the  £/0(£)  spectrum  and  exp 
[/</>(£)]  spectrum  for  each  component  are  separated  [Huang  et 
al. ,  1997].  In  the  following,  we  briefly  review  and  extend  the 
basic  ideas  of  EMD  method. 

It  is  noted  that  (5)  often  appears  as  a  WKB  solution  of  a 
wave  packet  to  a  wave  equation  that  includes  both  dispersive 
and  dissipative  phenomena  [e.g.,  Bender  and  Orszag,  1978] 

d2U 

-j^+HU=0,  (7) 


where  the  refractive  index  n  describes  wave-medium  interac¬ 
tion  and  it  can  be  complex.  We  define  the  complex  refraction 
index  by  n2  =  nf  +  inf.  In  general,  the  dispersive  phenom¬ 
enon  is  described  by  a  spatial  variation  of  nf,  whereas  the 
dissipation  effects  are  formulated  by  a  nonzero  nf  term.  When 
expressed  in  terms  of  (5),  the  WKB  solutions  of  (7)  are  given 
by  [e.g.,  Einaudi  and  Hines ,  1970;  Wait,  1981] 


tfo«)  = 


Cl,2 


'«»,(  t)2 
nM 


(8a) 


where  the  prime  denotes  the  derivative  (i.e.,  <f>'  =  d<t>/d f).  In 
this  paper,  we  are  mostly  interested  in  wind  fluctuations  as  a 
function  of  space.  Therefore  a  local  wavenumber  k(£)  =  <£'(£) 
for  the  complex  signal  w(f)  is  provided  by  (4).  When  expressed 
in  terms  of  amplitude  and  wavenumber,  the  complex  signal  (1) 
can  also  be  written  as 


U(&  =  U0(&  exp 


i  I  k(t)  dr 

Jo 


(5) 


An  optimal  choice  for  the  imaginary  function  iw(f)  that 
gives  an  unambiguous  local  wavenumber  is  the  Hilbert  trans¬ 
form  of  the  signal  u(£)  [Cohen,  1995],  defined  as 


i 

IT 


t 


m(t)  d T 


(6) 


where  the  integral  is  a  Cauchy  principle  value  and  H  is  the 
Hilbert  transform  operator.  The  complex  signal  thus  defined 
by  (1)  and  (6)  is  called  an  analytic  signal  belonging  to  «(f) 
[e.g.,  Bracewell,  1986;  Cohen,  1995].  Likewise,  t/0(f)  and  <£'(£) 
are  the  amplitude  and  wavenumber  of  the  analytic  signal  be¬ 
longing  to  the  original  signal  m(£).  Cohen  [1995]  shows  under 
certain  restricted  conditions  that  for  the  analytic  signal  defined 
by  (1)  and  (5)  there  is  no  overlap  between  the  wavenumber 
spectra  of  t/0(f)  and  exp  [*</>(£)]•  I n  other  words,  the  process 
of  constructing  an  analytic  signal  £/(£)  from  a  real  signal  m(£) 


k(£)  =  «,(€),  (8b) 

where  Cl  2  are  integration  constants  for  each  of  two  wave 
packets  propagating  in  opposite  directions.  The  posterior  con¬ 
dition  that  (5)  and  (8)  are  valid  solutions  of  (7)  is  a  slowly 
varying  envelope  of  the  wave  packet  for  each  wave  component 
[e.g.,  Einaudi  and  Hines,  1970;  Wait,  1981] 

1  d2U0\  1/2 

u,W I  <<K  =  nr  (9) 

If  we  define  the  left-hand  side  of  (9)  as  a  local  characteristic 
inverse  scale  K(£)  of  the  envelope  for  the  wave  packet,  then 
the  validity  condition  of  the  WKB  approximation  can  be  re¬ 
written  as 


1  d2(Jn\  1/2 

v.w\  -*r({)*K({)'  <10> 

Therefore,  in  order  that  an  analytic  signal  (5)  represents  a 
physically  meaningful  wave  packet  with  an  envelope  £/0(£)  and 
a  local  wavenumber  k(£),  the  local  characteristic  inverse  scale 
of  U0(£)  should  be  well  separated  from  the  local  wavenumber 
k(£)  of  the  analytic  signal.  Condition  (10)  can  be  considered  as 
an  acceptable  working  criterion  for  sufficient  separation  be¬ 
tween  the  £/(,(£)  spectrum  and  exp  [*<£(£)]  spectrum. 

Condition  (10)  requires  the  spectral  separation  at  local 
points  of  £  We  can  also  propose  the  following  global  condition 
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Figure  1.  Real  signals  of  a  superposition  of  two  sinusoidal  waves  defined  by  (12)  with  ax  =  1,  a  2  =  1.5, 
and  k1  =  1:  (top)  k2  =  0,9;  (middle)  k2  =  0,5;  and  (bottom)  k2  =  0.1.  The  dotted  lines  denote  0. 


for  spectral  separation  by  defining  a  wave  energy-weighted 
index  of  separation  x 


1  f  d 

l)-c 


Vo 

de 


Uld{ 


x  = 


4. 


T  KUodt 


«  1, 


(11) 


where  L  is  the  total  length  of  a  digital  signal. 

To  illustrate  the  above  conjecture  more  clearly  and  to  see 
how  the  EMD  method  generally  results  in  a  clear  separation 
between  K(€)  and  k(£)  for  each  empirical  mode,  we  now 
consider  a  well  known  example  of  a  superposition  of  two  si¬ 
nusoidal  signals  with  different  wavenumbers 

u(x)  =  ax  cos  {k^c)  +  a2  cos  ( k &).  (12) 


Without  loss  of  generality  we  first  fix  the  values  of  ax  =  1, 
a2  =  1.5,  and  k1  =  1.  Letting  k2  change  between  0  and  1,  we 
then  examine  U0(x),  K(x),  and  k(x).  In  Figure  1  we  show 
three  plots  of  u(x)  with  decreasing  values  of  k2  from  0.9  to  0.1. 
Although  the  high  frequency  oscillatory  feature  k1  =  1  exists 
in  all  three  signals,  Figure  1  shows  that  the  basic  character  of 
the  large-scale  variations  is  totally  different.  When  k2  **  kl9 
the  signal  resembles  a  wave  packet  with  slowly  varying  enve¬ 
lope,  while  when  k2  «  ky,  the  signal  represents  small  scale 
disturbances  riding  on  a  large-scale  background  field.  The  an¬ 
alytic  signal  belonging  to  (12)  can  be  derived  easily  based  on 
the  basic  formula  //[cos  (A:*)]  =  sin  ( kx )  [e.g.,  Bracewell , 
1986;  Karl,  1989].  The  local  amplitude  UQ(x)  and  the  local 
wavenumbers  k(x)  and  /£(*)  defined  by  (2),  (4),  and  (10)  for 
the  analytic  signal  are  [Cohen,  1995] 

U0(x)  =  [(a,  -  a2)2  +  4a1a2  cos2  {kjc)f2,  (13a) 


k(x)  =km-kd 


(13b) 


K(x)  = 


w 


U0(x)2 

\4ayaJt(a  y  —  a^)2cos(2kjc)  +  4a1a2cos4  (A;rf*r)]|1/2, 


(13c) 


where  km  and  k d  denote  the  average  and  the  difference  of  the 
two  basic  wavenumbers 


kd  = 


(14) 


Figure  2  shows  the  plots  of  U0(x),  k(x),  and  K(x)  for  the 
three  cases.  The  indices  of  separation  \  for  the  three  cases  are 
0.061,  0.432,  and  1.35,  respectively.  We  may  immediately  con¬ 
clude  the  following:  (1)  when  the  conditions  (10)  and  (11)  are 
well  satisfied  as  in  the  case  k2  =  0.9,  the  local  wavenumber 
k(x)  approximately  equals  the  fast  varying  basic  wavenumber 
kx  in  the  original  real  signal  and  varies  slightly  with*;  (2)  when 
(10)  and  (11)  are  not  satisfied,  the  local  wavenumber  k(x) 
departs  significantly  from  the  basic  wavenumber  k1  or  becomes 
negative.  To  see  how  the  index  of  separation  x  also  varies  with 
the  amplitudes  a  1  and  a2  of  basic  sinusoidal  waves,  we  show  in 
Figure  3  the  plot  of  Log10  (*)  as  a  function  of  two  parameters 
8  =  a2/a  i  and  e  =  kjkm ,  where 


#|e| 

X  IS"1  +  S  +  e(S-‘ -  S)| ’ 


(15) 


Figure  3  shows  that  the  index  of  separation  x  will  small 
when  the  amplitudes  of  the  two  basic  wave  components  differ 
significantly  (ax  <SC  a2  or  a 2  a±)  even  if  their  wavenum¬ 

bers  differ  significantly  (e  ~  1).  Clearly,  the  local  wavenumber 
of  the  analytic  signal  under  such  a  circumstance  is  approxi- 
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Figure  2.  The  amplitude  U0  (solid  lines)  and  local  wave  numbers  k(x)  (dashed  lines)  and  K(x)  (dash- 
dotted  lines)  of  the  analytic  signals  belonging  to  the  real  signals  of  Figure  1.  The  dotted  lines  denote  1. 


mately  equal  to  the  basic  wavenumber  that  has  a  dominant 
amplitude.  A  signal  with  a  {  a2ot  a2  corresponds  to 

a  major  wave  component  with  a  small  riding  wave  or  with  a 
small  asymmetric  distortion,  respectively. 

Returning  to  Figure  1  where  a  x  ~  a2,  we  find  that  the  fast 
varying  mode  corresponding  to  kx  ~  1  exists  visibly  in  all  the 
signals  no  matter  what  the  local  wavenumber  k(jc)  is.  The 
central  idea  of  the  EMD  method  is  to  separate  locally  the  fast 
varying  mode  from  the  background  slowly  varying  mode  by  a 
recursive  sifting  process  [Huang  et  aL,  1997].  The  sifting  pro¬ 
cess  extracts  a  signal  by  subtracting  the  mean  of  two  envelopes 
which  are  formed  by  smoothly  interconnecting  local  maxima 
and  local  minima  in  the  data  series,  respectively.  Each  such 
signal  is  called  an  Intrinsic  Mode  Function  (IMF).  Each  IMF 
is  subtracted  from  the  data  series,  and  the  process  is  repeated 
until  all  the  variability  in  the  data  series  is  decomposed  into  a 
number  of  these  IMFs.  Generally  speaking,  each  IMF  has  the 
form  of  a  wave  packet  which  satisfies  the  validity  condition 
(10).  For  a  more  precise  definition  of  the  IMF,  detailed  devel¬ 
opment  and  implementation  of  the  EMD  method,  Hilbert 
spectrum  analysis,  and  their  various  applications  in  analyzing 
nonlinear  and  nonstationary  time  series,  readers  are  referred 
to  the  work  of  Huang  et  al.  [1997].  When  applying  the  EMD 
method  to  the  present  example  (Figure  4),  we  recover  the  two 
basic  sinusoidal  components  as  their  IMFs  for  the  second  and 
third  cases  where  the  condition  (10)  and  (11)  are  not  well 
satisfied.  However,  the  IMF  for  the  first  case  is  the  signal  itself 
that  has  a  slight  local  wavenumber  modulation  and  a  significant 
amplitude  modulation  in  space. 

Since  the  derived  IMF  for  the  first  case  is  a  wave  packet  with 
an  amplitude  and  wavenumber  modulation,  then  the  EMD 
technique  may  be  particularly  useful  in  analyzing  and  inter¬ 
preting  dispersive-dissipative  wave-related  data.  We  know 
from  the  construction  that  the  signal  in  case  1  corresponds  to 
a  superposition  of  two  sinusoidal  waves  with  distinct  wavenum¬ 


bers  and  constant  amplitudes.  On  the  other  hand,  the  wave 
packet  IMF  derived  from  the  EMD  method  and  indicated  by 
(5)  and  (13)  suggests  that  there  exists  only  one  mode  in  the 
original  signal.  Although  the  two  expressions  are  identical 
mathematically,  they  may  have  different  physical  origins.  For 
example,  the  signal  in  the  first  case  can  be  produced  either  by 
two  wavemakers  that  have  distinct  frequencies  and  constant 
amplitudes  or  by  one  wavemaker  with  a  fixed  frequency  and  a 


Logl  O(delta) 

Figure  3.  Logarithm  of  index  of  separation,  Log10(*),  as  a 
function  of  Log10(S)  and  e. 
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Figure  4.  EMD  of  a  superposition  of  two  sinusoidal  waves:  IMFs  decomposed  from  real  signals  of  Figure 
1.  The  solid  lines  are  the  first  components  of  IMFs  corresponding  to  high  frequency  oscillations.  The  dashed 
lines  are  the  second  components  of  IMFs  corresponding  to  low  frequency  oscillations.  In  the  top  panel  there 
is  only  one  IMF  since  the  second  IMF  vanishes. 


slowly  variable  amplitude.  In  the  later  situation  the  medium  is 
inhomogeneous,  and  nr  is  slowly  vaiying  according  to  (13b). 

Since  both  expressions  of  the  signal  may  have  their  own 
physical  origins,  the  question  naturally  arises  as  to  which 
method  provides  a  better  spectral  decomposition.  Can  we  bet¬ 
ter  identify  the  actual  physics  behind  the  signal  by  using  the 
EMD  method  rather  than  the  traditional  Fourier  analysis?  The 
answer  is  generally  “yes”  for  nonstationary  wave-related  data 
since  IMFs  derived  from  the  EMD  method  provide  further 
information  on  spatial  variation  in  addition  to  the  amplitude 
and  wavenumber  for  each  mode.  Huang  et  at.  [1997]  present 
many  examples  extracted  from  numerical  results  of  nonlinear 
equation  systems  and  actual  data  to  show  the  power  of  the 
EMD  method.  In  the  next  section  we  define  a  local  power 
spectrum  and  its  marginal  distributions  which  can  be  applied  to 
analyzing  the  gravity  wave  characteristics  in  the  middle  atmo¬ 
sphere.  The  result  will  help  us  to  understand  what  mechanism 
is  mainly  responsible  for  the  observed  universal  spectrum  of 
gravity  waves  in  the  middle  atmosphere. 

Our  validity  condition  (10)  or  (11)  provides  a  precise  and 
quantitative  assessment  of  the  separability  between  the  U0({ ;) 
spectrum  and  exp  [*</>(£)]  spectrum.  It  also  provides  a  natural 
link  between  the  numerical  derived  IMFs  and  their  physical 
background  of  dispersive-dissipative  wave  packets.  Further¬ 
more,  such  a  condition  is  also  theoretically  crucial  for  calcu¬ 
lating  and  understanding  a  well-behaved  local  wavenumber  (or 
instantaneous  frequency).  Although  the  connection  between 
the  validity  condition  of  a  WKB  solution  for  a  wave  packet  and 
the  separability  condition  for  an  analytic  signal  is  natural  and 
inevitable,  previous  studies  in  signal  processing  have  appar¬ 
ently  missed  such  a  connection.  For  example,  Boashash  [1992] 
introduced  the  notion  of  monocomponent  signals  in  order  to 
find  a  well-behaved  instantaneous  frequency.  However,  his  in- 


vertibility  condition  excludes  most  of  the  signals  of  unmonoto- 
nous  property  which  do  have  a  well-defined  instantaneous  fre¬ 
quency.  We  have  noted  that,  historically,  a  WKB  solution  is 
mostly  used  for  solving  (7)  where  the  independent  variable  is 
space,  while  in  signal  processing  the  independent  variable  is 
usually  time.  To  conclude  this  section,  we  further  emphasize 
here  that  both  expressions  of  the  signal  for  the  case  1  in  Figure 
1  have  their  own  physical  backgrounds.  On  the  basis  of  the 
connection  between  the  separation  condition  of  the  spectra 
and  the  validity  condition  of  the  WKB  approximation  we  ex¬ 
pect  the  EMD  to  provide  a  natural  basis  for  the  analysis  of 
signals  which  contain  dispersive-dissipative  wave  phenomena. 

3.  Local  Power  Spectrum  and  Its 
Marginal  Distributions 

To  illustrate  how  the  EMD  method  may  be  helpful  to  un¬ 
derstanding  the  gravity  wave  characteristics  in  the  middle  at¬ 
mosphere,  we  need  to  introduce  a  distribution  function  that 
represents  the  energy  density  distribution  as  a  function  of 
space  and  wavenumber.  To  show  the  motivation  behind  its 
definition,  let  us  consider  an  idealized  wave  packet  of  velocity 
perturbation 

u(z ) 

exp  (z/20)  sin  [200  exp  (-z/21)]  6  £  z  ■<  50  km; 

exp  (2.5)  sin  [200  exp  (-z/21)]  50  s  z  <  65  ktn. 

(16) 

Figure  5  shows  the  velocity  profile  u(z),  its  Hilbert  transform 
H[U(z)]y  and  the  wave  packet  amplitude  of  tHe  analytic  signal 
U0(z).  Both  expression  (16)  and  Figure  5  suggest  that  the 


16,550 


ZHU  ET  AL.:  GRAVITY  WAVE  CHARACTERISTICS  IN  THE  MIDDLE  ATMOSPHERE 


-15.0  -10.0  -5.0  0.0  5.0  10.0  15.0 

Velocity  (m  s'1) 

Figure  5.  An  idealized  wave  packet  of  velocity  perturbation  w(z)  (solid  line),  its  Hilbert  transform  H[u  (z)] 
(dashed  line),  and  the  amplitude  of  the  analytic  signal  U0(z)  (dash-dotted  line). 


wave  packet  is  an  IMF  because  the  validity  condition  (10)  is 
well  satisfied.  Figure  6  shows  the  local  vertical  wavenumber 
m(z)  and  characteristic  inverse  scale  of  its  amplitude  envelope 
M(z)  of  the  analytic  signal  based  on  the  definitions  (4)  and 
(10),  respectively.  Two  spikes  in  m(z)  and  M(z)  at  50  km  are 
caused  by  the  discontinuity  of  du/dz  at  50  km.  Rapid  variations 
in  M(z)  close  to  boundaries  are  due  to  the  edge  effects  when 
performing  the  Fourier  transform  to  derive  H[u(z)]  [e.g., 
Karl,  1989].  Furthermore,  evaluation  of  the  second  derivative 
in  (10)  numerically  enhances  any  errors,  since  it  usually  in¬ 
volves  calculating  small  differences  of  large  terms.  In  practice, 
more  realistic  criteria  of  standard  deviation  and  indices  of 
orthogonality  [Huang  et  al ,  1997]  are  used  in  deriving  the 
iMFs. 

By  combining  Figure  5  and  Figure  6  one  can  establish  a 
one-to-one  correspondence  at  each  altitude  between  the  am¬ 
plitude  of  the  analytic  signal  U0(z)  (or  its  square  f/0(z)2)  and 
the  local  wavenumber  m(z).  The  distribution  of  energy  as  a 
function  of  wavenumber  and  space,  referred  to  as  the  Hilbert 
spectrum,  is  introduced  by  Huang  et  al.  [1997]  to  describe  the 
nonlinear  and  nonstationary  effects  of  various  signals.  In  the 
present  applications  of  analyzing  a  gravity  wave  spectrum  a 
more  relevant  quantity  is  the  energy  spectral  density  or  the 


made  between  the  EMD  method  and  other  time-frequency 
analysis  techniques  such  as  wavelet  analysis  and  the  Wigner 
distribution.  Wavelet  analysis  can  be  considered  as  an  im¬ 
proved  windowed  Fourier  transform  where  the  scales  are  ad¬ 
justable  but  the  basis  functions  are  specified  beforehand  [e.g., 
Daubec hies,  1992;  Kaiser ,  1994].  On  the  other  hand,  the  Wigner 
distribution  derives  internally  determined  basis  functions  from 
the  signal  itself,  but  the  transform  is  quadratic  rather  than 
linear  [Cohen,  1995].  The  EMD  method  possesses  both  merits 
of  linearity  and  internally  determined  basis  by  expressing  any 
signal  as  a  linear  superposition  of  a  finite  number  (/)  of  ana¬ 
lytic  signals  belonging  to  their  IMFs  [Huang  et  al.,  1997] 


u(z)  =  Re 


[■f 


2  U0](z)  exp  i  w/r)  dr 


(17) 


where  Re  {  }  denotes  the  real  part  of  the  complex  signal. 

Note  from  (17)  that  at  a  given  altitude  z  there  exists  only  a 
finite  numbers  of  wavenumbers  Wy(z).  Therefore  it  is  conve¬ 
nient  to  discretize  both  space  and  wavenumber  variables  ac¬ 
cording  to 


power  spectrum,  F(m),  which  often  shows  a  near  —3  power 
law  at  large  m .  In  this  section  we  first  introduce  a  distribution 

zk  —  Zo  +  (k  —  1/2)A  z 

k  =  \,  2,  •  • 

• ,K . 

(18a) 

function  which  we  call  the  local  power  spectrum.  We  then 

m(+i  =  mt  +  Am, 

€  =  1,2,  •• 

' ,  L. 

(18b) 

suggest  that  different  physical  mechanisms  for  the  gravity  wave 
spectrum  may  lead  to  different  observed  structure  in  the  local 
power  spectrum  even  though  the  corresponding  standard 
power  spectra  F(m)  are  indistinguishable  for  different  pro¬ 
cesses. 

In  the  work  of  Huang  et  al.  [1997],  comparisons  are  also 


In  the  computational  algorithm  the  altitude  z  will  be  assigned 
toz*  if  zk  —  Az/2  <  z  <  zk  +  Az/2.  A  similar  procedure  is 
also  applied  to  discretizing  the  wavenumber  m}  according  to 
the  limits  in  (18b).  The  local  power  spectrum  at  a  space- 
wavenumber  point  {zk,  m€)  is  defined  as 
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Local  Vertical  Wavenumber  (km1) 

Figure  6.  Local  vertical  wavenumbers  m(z)  (solid  line)  and  M{z)  Of  the  analytic  signal  in  Figure  5  based 
on  the  definitions  (4)  and  (10),  respectively. 


1 

P(zh  mt)  =  j- 


2 

;=i 


[£/0j(zfc;  /w/z^))]2 
(A  m() 


8{m(,  m). 


(19)  G(zk)  =  2k*Z  P{zk,  me)Amt  k  =  1,  2,  •  •  • ,  K.  (22b) 


where  8(me,  m})  is  the  Kronecker  6  function 


8(m€,  m})  = 


m(  =  m}{zk ); 
m<  #  m}{zk ). 


(20) 


In  (19)  we  have  written  U0j(zk)  as  U0j(zk;  mj(zk))  to  indicate 
that  the  amplitude  is  evaluated  at  z  —  zk  where  the  corre¬ 
sponding  vertical  wavenumber  is  m;-.  Clearly,  the  local  power 
spectrum  P(zk,  m€)  is  nonnegative  and  has  the  conventional 
units  of  a  spectral  density,  in  this  case  variance  per  unit  wave- 
number.  The  factor  1/L  in  the  definition  makes  P(zk ,  m€) 
nearly  constant  for  different  resolutions  of  the  wavenumber 
(e.g.,  L*  =  2L  if  Am*  =  Am€/2).  P(zk,  m€)  can  be  con¬ 
sidered  as  a  two-dimensional  discrete  probability  density  func¬ 
tion  with  the  normalization  constant 


N*  =  2K  ^  ^  P^Zk’ 


(21) 


Here  the  marginal  distribution  F{m€)  is  a  measure  of  power 
spectrum  at  the  wavenumber  while  G(zk)  is  the  variance 
at  the  altitude  zk.  Under  the  condition  of  constant  U0j  and  my 
such  as  in  the  last  two  cases  of  Figure  1  the  normalization 
constant  is  the  average  energy  of  the  real  signal  w(z);  that  is, 


1  * 

Ntt  =  £  2  u(zk)2, 


(23) 


and  the  marginal  distribution  F(me)  is  the  power  spectrum  of 
u(z)  derived  from  the  Fourier  transform.  In  Plate  1  and  Figure 
7  we  show  the  local  power  spectrum  P{zk ,  m€)  and  its  mar¬ 
ginal  distribution  F(m€),  respectively,  for  the  wave  packet 
signal  (16)  where  it  has  only  one  IMF  (/  =  1).  To  demon¬ 
strate  that  the  EMD  method  is  superior  to  the  wavelet  analysis 
in  this  case,  we  also  present  in  Plate  1  the  local  power  spectrum 
of  wavelet  analysis  by  replacing  the  amplitude  U0  in  (19)  with 
the  continuous  wavelet  transform  [Daubechies,  1992] 


The  two  marginal  distributions  associated  with  P(zk,  m€)  can 
be  defined  as 

a: 

F(mt)  =  2  P(zki  m()  €  =  1,  2,  •  •  • ,  L;  (22a) 

k=  1 


C0(zk,  mt)  =  -j=  T  u(r)i/r*(-  dr,  (24) 

where  a  =  'irim€  is  the  scale  parameter  and  the  asterisk 
denotes  complex  conjugation.  The  standard  Morlet  wavelet 
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Figure  7.  Marginal  distribution  F(m€)  derived  from  EMD  method  (solid  line)  and  power  spectrum  derived 
from  the  Fourier  transform  (dashed  line)  for  the  idealized  wave  packet  of  Figure  5. 


iKt)  =  Ce-’^V’"'  -  e-^14)  with  C  =  V2/a  and  a  =  4 
has  been  used  in  Plate  1.  We  can  see  from  Plate  1  that  although 
the  wavelet  analysis  provides  space-wavenumber  information 
for  the  signal,  the  spread  of  the  actual  local  wavenumber  to 
neighboring  heights  is  unavoidable  because  C0(zk,  me)  is  a 
continuous  function  of  space  and  scale.  On  the  basis  of  the 
construction  of  our  idealized  waveform  (16)  we  can  conclude 
that  such  a  spread  represents  an  energy  leakage  in  the  wave- 
number  domain  due  to  a  less-than-ideal  choice  for  the  basis 
function.  Mathematically,  this  is  also  due  to  the  admissibility 
condition  of  the  mother  wavelet  ^(r)  that  yields  a  norisingular 
property  of  i//(t)  at  r  =  0.  More  examples  of  comparison 
between  the  EMD  method  and  the  wavelet  analysis  are  pre¬ 
sented  by  Huang  et  al.  [1997]. 

From  the  definition  of  the  local  power  spectrum  we  know 
that  each  IMF  corresponds  to  one  curve  in  the  two- 
dimensional  contour  plotP(z*,  me).  As  the  number  /  of  IMFs 
increases,  one  expects  a  smoother  contour  plot  P{zk,  m€)  that 
may  provide  a  continuous  and  positive  distribution  function. 
However,  in  real  applications,  /  usually  is  a  small  number.  For 
example,  to  test  our  sifting  program,  we  applied  the  EMD 
method  to  a  data  series  believed  to  be  most  irregular  and 
unpredictable:  a  30-year  long  record  of  the  composite  index  of 
The  New  York  Stock  Exchange  with  7592  trading  days.  Apply¬ 
ing  the  EMD  method  to  the  time  series  yields  only  10  IMFs 
with  timescales  ranging  from  daily  variability  to  an  interdec- 
adal  oscillation,  Therefore,  in  principle,  P(zk,  me)  is  discrete 
in  wavenumber  (me). 

It  should  be  noted,  however,  that  a  signal  or  a  data  series  can 


be  considered  as  a  realization  (or  a  sample)  of  a  physical 
system.  In  many  applications  the  actual  distribution  function 
P(zk,  m€)  of  a  physical  system  can  be  derived  by  an  ensemble 
average  of  many  independent  realizations.  Under  these  cir¬ 
cumstances,  P(zkim€)  may  show  its  continuous  nature  in  both 
space  and  wavenumber  if  the  actual  physical  system  behaves  in 
such  a  manner.  In  the  present  application  to  the  vertical  wave- 
number  spectrum  of  gravity  waves  in  the  middle  atmosphere, 
P(zk1  m£)  is  expected  (according  to  some  models)  to  be  con¬ 
tinuous  due  to  sporadic  wave  generation,  dispersive  wave  prop¬ 
agation,  and  nonlinear  wave-wave  interactions.  Other  models, 
however,  suggest  it  may  be  somewhat  more  discrete  in  certain 
circumstances  [e.g.,  Alexander ,  1996;  Warner  and  McIntyre , 
1996]. 

When  there  is  a  frequency  or  an  amplitude  modulation  in  an 
IMF,  the  marginal  distribution  FEMD(m)  and  the  Fourier 
spectrum  FFoiirier(m)  will  be  different.  For  example,  from  (13) 
we  find  FEMD(m)  for  the  IMF  in  Figure  1  to  be 

a\  -  af 

iWK)=-^— p,  (25) 

where  the  wavenumber  k  varies  continuously  within  a  wave- 
number  domain  of  a  finite  support  [km  +  kd(a2  -  ax), 
km  +  kd(a2  +  ax)\.  Outside  such  a  finite  wavenumber  sup¬ 
port,  FBMD(m)  vanishes.  On  the  other  hand,  the  Fourier  spec¬ 
trum  FFourier(m)  of  the  same  signal  consists  of  two  spectral 
lines  located  at  k  =  k1  and  k  =  k2 ,  respectively.  As  we  have 
already  discussed  in  the  previous  section,  either  form  of 
Femd(k)  or  FFourier(m)  may  correctly  correspond  its  own 
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physical  process.  One  usually  cannot  identify  the  actual  phys¬ 
ical  process  based  on  differences  between  FEMD(m)  and 
irFourier(m)  alone.  This  example  also  shows  the  importance  of 
choosing  different  spectral  analysis  tools  in  analyzing  the  data. 
The  Fourier  spectrum  resolves  two  harmonics  since  it  assumes 
that  the  series  is  a  superposition  of  infinite  constant  amplitude 
harmonics.  Thus  the  only  way  for  it  to  represent  a  slowly 
varying  amplitude  is  to  split  it  up  into  two  harmonics  which 
beat  when  retransformed.  The  EMD  method  avoids  this  and 
gives  a  direct  picture  of  amplitude  modulation  of  the  wave 
packet.  However,  conversely,  the  EMD  method  cannot  resolve 
two  closely  spaced  harmonics  due  to  the  acceptable  separabil¬ 
ity  criterion  (10)  or  (11).  In  the  real  applications,  such  as  in  the 
present  situation,  when  a  data  series  is  associated  with  physical 
processes  which  can  be  approximated  by  the  wave  equation  (7), 
the  EMD  method  is  expected  to  yield  a  better  decomposition. 

For  a  real  signal  with  multi-IMFs  and  a  finite  length  its 
FEMD(m)  will  still  be  a  continuous  function  of  a  finite  wave- 
number  support  because  the  local  wavenumbers  are  'continu¬ 
ous  functions  of  space.  On  the  other  hand,  the  Fourier  spec¬ 
trum  of  a  finite  length  signal  FFourier(m)  will  be  unlimited  in 
wavenumber  space.  The  cutoff  wavenumbers  are  determined 
mainly  by  how  the  continuous  signal  is  sampled  [e.g.,  Karl, 
1989;  Harris ,  1978].  It  is  a  notoriously  difficult  problem  to  find 
a  positive  space -wavenumber  (or  time-frequency)  distribution 
that  satisfies  the  specified  marginal  distributions  [Cohen,  1995]. 
Because  our  marginal  distributions  (22)  are  defined  by  the 
nonnegative  local  power  spectrum,  the  required  condition  of  a 
two-dimensional  probability  density  is  satisfied  automatically. 
Although  our  marginal  distributions,  such  as  FBMD(m),  are 
different  from  the  marginal  distributions  specified  in  the  usual 
time-frequency  analysis  [Cohen,  1995],  such  as  F¥ouTtGT(m), 
the  marginal  distributions  are  physically  meaningful  as  we  dis¬ 
cussed  in  the  last  section.  Therefore  the  local  power  spectrum 
P{zk,  m€)  is  clearly  an  appropriate  distribution  function. 

It  is  also  noted  that  the  Fourier  spectrum  FFourier(m)  de¬ 
rived  in  turbulence  studies  also  represents  an  ensemble  aver¬ 
age  of  many  independent  realizations.  However,  for  homoge¬ 
neous  turbulence  where  the  flows  are  considered  to  be 
stationary  in  time  the  ergodic  theorem  ensures  that  the  ensem¬ 
ble  average  is  the  same  as  the  average  over  a  long  data  series 
[e.g.,  Panchev ,  1971;  Frisch,  1995].  Examining  Figure  5  clearly 
shows  that  the  wave  packet  is  not  a  stationary  signal.  Therefore 
additional  cautions  should  be  exercised  to  interpret  FFourier(m) 
when  it  is  derived  from  only  a  few  data  series. 

4.  Gravity  Wave  Characteristics 
in  the  Middle  Atmosphere 

On  returning  to  Figure  7  we  are  now  in  a  position  to  dem¬ 
onstrate  some  of  the  usefulness  of  the  EMD  method  and  the 
local  spectrum  P(zk,  m€)  in  the  analysis  of  gravity-wave  data 
in  the  middle  atmosphere.  We  have  already  discussed  the  dif¬ 
ferences  between  FBMD(m)  and  FFourier(m)  in  terms  of  then- 
physical  interpretations.  The  most  significant  difference  is  that 
they  use  a  different  set  of  basis  functions  to  decompose  the 
data,  as  shown  in  Figure  7.  However,  within  the  major  part  of 
their  common  wavenumber  domain,  FBMD(m)  is  very  close  to 
FFounGT(m)  since  both  functions  represent  the  wave  energy  per 
wavenumber.  In  the  present  example  of  the  idealized  wave 
packet  shown  in  Figure  5,  both  FBMT>(m)  and  F¥ouTieT(m) 
exhibit  a  —3  slope  in  their  Log  (F)-log  ( m )  plots.  This  is 
consistent  with  the  vertical  wavenumber  spectra  derived  from 


many  observations  [e.g.,  VanZandt ,  1982;  Tsuda  et  al.,  1989; 
Wu  and  Widdel,  1989, 1992;  Senftand  Gardner,  1991].  Although 
the  -3  slope  of  FEMD(m)  and  FFourier(m)  in  the  present 
example  can  be  considered  as  coincidental,  the  idealized  wave 
packet  shown  in  Figure  5  resembles  observations  in  several 
respects.  First,  observations  also  show  that  the  amplitude  of 
internal  gravity  waves  increases  exponentially  with  altitude 
[e.g.,  Hirota  and  Niki,  1985]  but  that  wave  amplitudes  tend  to 
be  constant  with  height  in  the  mesosphere  [see,  e.g.,  Tsuda  et 
al.,  1994,  Figure  7].  The  first  observation  is  consistent  with  the 
notion  that  the  density  effect  of  the  background  atmosphere 
amplifies  the  upward  propagating  gravity  waves  [Hines,  1960], 
whereas  the  second  observation  is  consistent  with  the  per¬ 
ceived  “saturation”  of  wave  amplitude  at  mesosphere  heights. 
Second,  changes  of  refraction  index  associated  with  back¬ 
ground  wind  shear  and  varying  static  stability  also  alters  the 
vertical  wavenumber  [e.g.,  Bleistein,  1984;  Bretherton,  1966; 
Zhu,  1987;  Marks  and  Eckermann,  1995]. 

The  near  -3  slope  of  the  vertical  wavenumber  spectrum  for 
the  gravity  waves  in  the  middle  atmosphere  was  first  noted  by 
VanZandt  [1982].  He  suggested  that  such  a  universal  spectrum 
in  the  atmosphere  might  be  similar  to  the  universal  spectrum 
introduced  by  Garrett  and  Munk  [1972,  1975]  to  describe  oce¬ 
anic  internal  waves.  Since  then,  many  different  theories  have 
been  proposed  to  explain  the  apparent  universal  -3  power  law 
for  the  vertical  wavenumber  spectrum  [e.g.,  Dewan  and  Good, 
1986;  Smith  etai,  1987;  Weinstock,  1990;  Hines,  1991;  Gardner , 
1994;  Zhu,  1994].  We  can  divide  different  mechanisms  into  two 
arbitrary  categories:  linear  versus  nonlinear.  A  linear  mecha¬ 
nism  usually  suggests  that  there  is  a  characteristic  (or  domi¬ 
nant)  wave  component  with  the  amplitude  growing  with  alti¬ 
tude  due  to  the  density  effect.  However,  the  characteristic  wave 
component  grows  at  a  rate  smaller  than  the  inverse  of  the  air 
density  due  to  the  existence  of  damping  mechanisms  such  as 
wave  breaking  and  radiative  damping.  An  ensemble  average  of 
all  characteristic  wave  components  forms  a  near  —3  power 
spectrum.  The  amplitude  of  a  characteristic  wave  is  also 
strongly  controlled  by  the  different  generation  mechanisms 
[e.g.,  Lindzen,  1984;  Zhu  and  Holton,  1987].  In  nonlinear  mech¬ 
anisms  all  wave  components  interact  in  a  way  that  in  some 
models  cascades  the  larger-scale  waves  to  smaller  scales  [e.g., 
Dewan,  1994]. 

Even  though  both  marginal  distributions  of  linear  and  non¬ 
linear  mechanisms  may  yield  the  same  form  of  power  law 
fEMD(m)  ~  m~p  with  p  ~  3,  its  local  power  spectra 
P(zk9  m€)  will  exhibit  different  structures  in  terms  of  their 
spatial  and  wavenumber  distributions.  In  the  case  that  a  single 
wave  component  dominates  the  whole  spectrum  one  expects 
P{zk,  m€)  to  be  sharply  peaked  along  one  line  as  in  Plate  1. 
Furthermore,  the  strengths  and  mechanisms  that  excite  gravity 
waves  at  different  locations  and  altitudes  differ  significantly.  As 
a  consequence,  at  a  given  altitude  zk,  P(zk,  m€)  may  not  have 
an  m~p  power  spectrum  as  a  function  of  m€.  On  the  other 
hand,  when  the  nonlinear  cascade  process  becomes  dominant, 
one  may  expect  a  universal  m  ~p  power  spectrum  inm^  that  is 
independent  of  location  and  altitude. 

To  test  which  mechanism  is  closer  to  the  observations,  we 
present  in  Plate  2a  the  local  power  spectra  P(zk,  m€)  calcu¬ 
lated  from  38  horizontal  wind  profiles  (u,  v)  obtained  during  the 
Dynamics  Adapted  Network  for  the  Atmosphere  (DYANA) 
campaign,  carried  out  in  January  through  March  1990  in  the 
northern  hemisphere  [Offermann,  1994].  All  38  profiles  were 
derived  by  tracking  passive  falling  spheres  released  from  rock- 


LoglO(m) 

Plate  1.  Logarithm  of  local  power  spectrum,  Log10  [P(zk,  me)],  of  the  idealized  wave  packet  of  Figure  5: 
(top)  P(zk,  m()  is  calculated  by  (19)  from  the  EMD  method;  (bottom)  P(zk,  m()  is  calculated  by  (19)  from 
the  wavelet  analysis  with  U0  in  (19)  replaced  with  C0  in  (24).  The  color  scheme  is  made  to  “saturated”  below 
the  value  of  -5  to  accentuate  significant  wave  energy  in  the  neighborhood  of  the  local  wavenumber. 
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Plate  2.  Logarithm  of  local  power  spectrum,  Log10  [P(zk,  m()],  of  the  observed  horizontal  wind  fields: 
(top)  from  DYANA  campaign  (69°N,  16°E  and  44°N,  1°W)  and  (bottom)  from  long-term  rocket  soundings  at 
Thule  (77°N,  69°W). 


ZHU  ET  AL.:  GRAVITY  WAVE  CHARACTERISTICS  IN  THE  MIDDLE  ATMOSPHERE 


16,555 


ets  at  an  apogee  —95  km,  with  data  available  in  an  altitude 
region  between  33  and  91  km.  Fifteen  rockets  were  launched  at 
Biscarosse  (France;  44°N,  1°W),  and  23  rockets  were  launched 
at  Andpya  (Norway;  69°N,  16°E).  It  is  found  from  Plate  2a  that 
the  local  power  spectra  P(zk,  m€)  is  peaked  at  a  wavenumber 
m  *  which  varies  with  altitude 


m*(z) 


2.0  exp  [-(z  -  40)/19.5]  z  <  78  km; 

0.28,  z  >  78  km. 


(26) 


Here  we  define  m  *  as  the  characteristic  vertical  wavenumber 
at  which  the  local  power  spectra  reach  their  maxima.  Such  a 
definition  is  also  consistent  with  the  m  *  defined  by  Smith  et  al. 
[1987]  where  m*  is  also  the  lower  cutoff  wavenumber  for  a 
m~ 3  region  in  the  power  spectrum  of  gravity-wave-induced 
horizontal  velocity  perturbations. 

A  striking  feature  shown  in  Plate  2a  is  a  disjointed  secondary 
maximum  in  P(zk,  me)  above  60  km  with  its  peak  vertical 
wavenumber  mf  1  order  of  magnitude  greater  than  the  pri¬ 
mary  m  *  defined  by  (26).  The  clear  separation  between  m  * 
and  mf  shown  in  P(zk,  m€)  suggests  that  there  exists  two 
different  physical  processes  above  60  km.  Because  of  the  at¬ 
mospheric  density  effect,  the  response  function  of  a  falling 
sphere  to  high  wavenumber  wind  fluctuations  decreases  rapidly 
with  the  altitude  [Hass  and  Meyer ,  1987].  Therefore  the  sec¬ 
ondary  maximum  in P(zk,  m€)  above  60  km  will  correspond  to 
very  strong  gravity  wave  activities  if  it  were  due  to  wind 
fluctuations.  We  believe  that  this  secondary  maximum  in 
P(zk ,  m€)  is  due  to  sphere  vibrations  caused  by  shock  waves 
emitted  from  the  falling  spheres  [e.g.,  Lighthill ,  1978;  Landau 
and  Lifshitz ,  1987]. 

It  can  be  easily  shown  that  a  free-falling  sphere  reaches  the 
speed  of  Mach  1  (speed  of  sound  **300  ms-1  around  the 
winter  mesopause)  at  a  traveling  distance  of  yHI  2  when  ini¬ 
tially  released  at  rest,  where  y  =  cplc  v  =  1 .4  is  the  ratio  of  the 
specific  heats  and  H  is  the  scale  height.  At  a  certain  stage  the 
falling  sphere  will  reach  a  terminal  velocity  that  decreases  with 
the  decreasing  altitude  because  the  air  drag  increases  [Miller, 
1969].  Eckermann  and  Vincent  [1989]  estimated  that  the  ter¬ 
minal  velocity  at  60  km  is  about  Mach  0.5  (150  ms'1).  There¬ 
fore  we  expect  spheres  to  fall  supersonically  above  —65  km. 
For  a  sphere  falling  at  a  speed  greater  than  Mach  1  the  shock 
waves  are  emitted  downward  within  a  cone  with  an  acute  angle 
of  cos-1  (1/M),  where  M  is  the  Mach  number.  The  emission 
of  shock  waves  may  cause  the  sphere  to  experience  a  jolt  at 
high  frequency  which  results  in  a  high  wavenumber  oscillation 
in  the  retrieved  vertical  profile  at  that  altitude.  Between  60  and 
65  km  the  sphere  could  still  suffer  a  forced  oscillation  due  the 
catch-up  of  the  shockwaves  previously  emitted  above  65  km.  It 
should  be  noted  that  because  of  the  atmospheric  density  effect, 
the  secondary  peak  wavenumber  m1"  tilts  with  height  in  a  form 
similar  to  m*.  Such  an  unrealistic  feature  about  the  wind 
profiles  will  be  difficult  to  identify  by  conventional  Fourier 
analysis. 

To  quantitatively  examine  which  mechanism  is  more  respon¬ 
sible  for  a  possibly  observed  near  -3  slope  in  the  Fourier 
analysis  [e.g.,  Zhu,  1994],  we  show  in  Figure  8  the  marginal 
distributions  at  different  altitude  ranges.  We  see  from  the  fig¬ 
ure  that  the  marginal  distribution  for  the  entire  region  from  35 
to  85  km,  shows  an  overall  near  -3  slope  except 

around  4  km'1  where  the  localized  energy  density  bump  due  to 
the  shock  waves  produces  to  a  -2  slope.  Again,  without  the 
knowledge  of  the  local  power  spectrum  it  will  be  difficult  to 


realize  two  distinct  processes  that  are  responsible  for  the  ap¬ 
parently  smooth  —2  slope  in  the  marginal  distribution  or  the 
power  spectrum  from  Fourier  analysis.  By  looking  into  indi¬ 
vidual  contributions  from  the  local  power  spectra  to  the  mar¬ 
ginal  distribution  we  can  identify  some  individual  physical  pro¬ 
cesses  that  are  responsible  for  the  form  of  the  marginal 
distribution,  F^m). 

Next,  we  calculate  the  marginal  distribution  for  the  entire 
region  but  contributed  spectrally  only  from  the  following  wave- 
number  band,  Fm*(m), 

Log  10(m*(z))  -  0.25  <  Log10(m)  <  Logi0 (m*(z))  +  0.25, 

(27) 

where  ra*(z)  is  defined  by  (26).  The  result  is  also  plotted  in 
Figure  8.  We  see  that  when  the  vertical  wavenumber  m  <  0.1 
km-1  (or  vertical  wavelength  Lz  >  63  km),  the  two  marginal 
distributions  are  almost  identical  indicating  that  characteristic 
gravity  waves  parameterized  by  a  single  wavenumber  contain 
most  of  the  energy.  Such  a  picture  is  consistent  with  the  “quasi- 
monochromatic”  parameterization  of  gravity  wave  processes 
advocated  by  Yamanaka  and  Fukao  [1994]  where  the  dominant 
vertical  wavenumbers  have  form  similar  to  (26).  Furthermore, 
the  marginal  distribution  Fm*(m)  is  much  closer  to  an  m~3 
power  law  than  F^m).  Returning  to  Plate  2a,  we  now  can 
clearly  identify  that  the  small  bump  of  a  local  m  ~2  power  law 
in  Fall(m)  is  caused  by  a  secondary  maximum  in  P{zk,  m€) 
above  60  km  that  has  a  vertical  wavenumber  raf  ~  2  km'1 
(Lz  —  3  km).  To  exclude  the  nonstationary  effects  of  the 
atmospheric  density  and  stratification,  we  calculate  and  plot  in 
Figure  8  the  marginal  distributions  of  3  km  slabs  centered  at 
different  altitudes  Fslah(m).  After  excluding  the  unphysical 
bumps  above  60  km,  we  see  that  when  m  >  3  km'1  (Lz  <  2 
km),  the  marginal  distributions  for  the  slabs  generally  exhibit  a 
power  law  distribution  m  ~p  with p  -  2  —  4.  This  result  is  also 
consistent  with  the  wavelet  spectrum  analysis  of  the  lower 
stratospheric  temperature  and  wind  observations  [Sato  and 
Yamada,  1994].  When  m  <  3  km'1  ( Lz  >  2  km),  the  slab 
marginal  distributions  show  a  distinct  feature  of  systematic 
shifting  of  the  peak  Fslab(m )  toward  smaller  vertical  wavenum¬ 
ber  as  the  altitude  increases.  The  shifting  is  consistent  with 
variation  of  the  characteristic  vertical  wavenumber  along  the 
altitude,  as  indicated  in  Figure  8.  However,  there  appears  no 
universal  power  law  that  can  fit  Fslah(m)  for  m*  <  m  <  3 
km'1  even  though  a  -3  power  law  exists  in  Fall(wi)  in  this 
wavenumber  range. 

In  Plate  2a  we  also  show  the  local  power  spectra  P(zk,  me) 
calculated  from  175  horizontal  wind  profiles  (w,  v)  from  me¬ 
teorological  rocket  soundings  at  Thule  (77°N,  69°W).  The  data 
cover  the  altitude  range  from  20  to  60  km  with  1  km  vertical 
resolution  [Eckermann  et  al.,  1994].  Here,  we  see  again  a  rel¬ 
atively  concentrated  P(zk,  m€)  around  m  =  1  km'1  (Lz  = 
6.3  km).  Furthermore,  since  more  wind  profiles  have  been 
used,  P(zk,  m€)  is  smoother  than  in  Plate  2a.  However,  be¬ 
cause  of  the  lower  resolution  of  the  data,  many  features  of 
P(zk,  m€)  at  higher  vertical  wavenumbers  may  have  been  lost. 

5.  Hodographic  Analysis  of  Characteristic 
Gravity  Waves 

We  have  established  in  section  2  a  close  connection  between 
the  WKB  wave  packet  solutions  of  a  dispersive -dissipative 
physical  system  and  IMFs  derived  from  the  EMD  method.  In 
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Figure  8.  Marginal  distributions  at  different  altitudes  calculated  from  Plate  2a.  Two  solid  lines  are  the 
marginal  distributions  in  the  entire  region  ranging  from  35  to  85  km.  The  one  with  the  open  squares  is 
calculated  by  setting  P(zk ,  m€)  =  0  outside  the  spectral  band  of  characteristic  vertical  wavenumber  as 
defined  by  (27).  The  dashed  lines  are  the  marginal  distributions  within  a  3  km  slab  centered  at  the  altitude  as 
indicated.  The  thin  solid  straight  lines  represent  either  m  “3  or  m  “2  power  law.  Five  vertical  parallel  bars  mark 
the  characteristic  vertical  wavenumbers  m  *  at  different  altitudes.  All  the  curves  except  the  one  at  40  km  have 
been  consecutively  moved  up  by  3  orders  of  magnitude  to  avoid  clustering. 


the  work  of  Huang  et  al.  [1997],  several  examples  are  presented 
to  show  that  many  derived  IMFs  have  their  own  physical  in¬ 
terpretations.  Here  we  show  for  the  DYANA  campaign  data 
that  one  IMF  can  be  identified  as  a  major  component  of 
upward  propagating  inertia  gravity  waves  in  the  middle  atmo¬ 
sphere. 

In  falling  sphere  experiments  the  wind  and  temperature 
profiles  are  extracted  from  the  retrieved  atmospheric  density 
profiles.  To  exclude  background  mean  and  planetary  wave 
scale  perturbations  and  to  account  for  the  strong  nonstation- 
arity  of  the  instrument  response  function  [Hass  and  Meyer , 
1987],  we  first  employ  a  band-pass  filter  to  decompose  a  wind 
profile,  u  or  v ,  into  background  and  perturbation  components 

u(z)  =  U(z)  +  w'(z).  (28) 


u'(z)  =  2  11,(2) 

2  Uo,(z)  exp  fij  m,( t) 


=  Re 


dr 


(29) 


whereas  for  the  background  component  U(z ),  EMD  method 
is  used  to  derive  one  IMF  denoted  as  ua 


U{z)  =  ub(z)  +  ug(z)  =  ub(z) 


+  Re 


U0g(z)  exp 


mg(r)  dr 


(30) 


We  use  a  101  points  Savitzky-Golay  filter  corresponding  to  a  2 
km  window  width  [Press  etal. ,  1992,  p.  644].  The  EMD  method 
is  applied  to  the  perturbation  uf(z)  to  derive  all  of  its  IMFs 


Here  the  residual  in  the  low  wavenumber  component  ub(z)  is 
considered  as  the  background  mean  and  planetary  wave  scale 
perturbations.  We  apply  the  same  procedure  to  the  meridional 
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figure  9.  (a)  Zonal  wind  u  observed  at  Biscarosse  (France;  44°N,  1°W)  and  its  decomposition  by  (28)  and 
(29):  ub,  iig,  and  w\  In  the  last  strip,  solid  line  and  dashed  line  denote  vg  and  H(ug),  respectively,  (b)  Five 
IMFs  for  the  high-pass  zonal  wind  u '  shown  in  Figure  9a. 


components  v(z).  Our  local  power  spectra  presented  in  Plate  and  amplitude  modulations.  It  is  well  known  that  at  the  mid 
2  are  calculated  based  on  all  the  IMFs  ug(z)  and  2  Uj(z).  and  high  latitudes,  horizontal  wind  perturbations  of  a  pure 

In  Figures  9  and  10  we  present  two  examples  of  decompos-  gravity  wave  component  show  elliptical  rotation  with  height 

ing  the  original  wind  profiles  into  various  components  ub,  ug ,  [e.g.,  Hirota  and  Niki ,  1985].  This  is  because  the  polarization 

ur ,  and  Uj.  We  can  see  from  figures  that  ug  contains  most  of  relation  of  inertia  gravity  waves  renders  a  constant  phase  shift 

the  wave  energy  for  the  disturbances  of  scales  less  than  those  between  ug  and  vg  components  [e.g.,  Gossard  and  Hooke , 

of  ub.  The  structures  of  ug  in  these  two  examples  resemble  1975].  Furthermore,  the  Hilbert  transform  of  a  signal  also 

upward  propagating  gravity  waves  with  certain  wavenumber  produces  a  90°  phase  shift  from  the  original  data  [e.g., 
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Figure  10.  Same  as  Figure  9,  except  for  a  wind  profile  at  Andoya  (Norway;  69°N,  16°E). 


Bracewell,  1986;  Karl,  1989].  In  the  last  strips  of  Figures  9a  and  mann  and  Vincent  [1989]  suggested  that  the  Stokes  parameters 

10a  we  also  show  vg  an6H(ug)  profiles.  Clearly,  vg  and  H(ug)  could  be  calculated  spectrally.  Recently,  Eckermann  [1996] 

exhibit  a  near-constant  phase  difference  in  the  majority  of  the  investigated  the  relationships  among  the  Stokes  parameters, 

region,  strongly  suggesting  that  ( ug ,  vg)  are  representing  a  rotary  spectra,  and  cross-spectral  methods.  Given  the  Stokes 

near  monochromatic  inertia  gravity  wave.  parameters,  some  gravity  wave  characteristics  such  as  the 

A  convenient  and  quantitative  way  to  describe  the  wave  phase  difference  8g  and  the  major  axis  orientation  rg  can  be 

propagation  characteristics  is  to  use  four  Stokes  parameters  derived  [Eckermann  and  Vincent,  1989] 

(/,  D,  P,  Q).  The  Stokes  parameters  were  first  introduced  for 

studying  the  gravity  waves  by  Vincent  and  Fritts  [1987].  Ecker-  dg  -  tan-1  ( QlP ), 


(31) 
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Figure  11.  Correlation  coefficient  Corr  [ug,  H(  vg)\  and  the  major  axis  orientation  rg  derived  for  all  38  wind 
profiles. 


t,-  1/2  tan'1  (P/D).  (32) 

In  the  above  the  axial  anisotropy  D  and  the  linear  polarization 
P  can  be  straightforwardly  calculated  from  velocity  profiles 
[j Eckermann ,  1996] 

D  =  -  7ff,  (33) 

P  =  2UgVg ,  (34) 

where  overbars  denote  height  averaging  of  the  velocity  profiles 
over  an  altitude  range  longer  than  local  vertical  wavelengths. 
With  the  help  of  the  Hilbert  transforms  of  ug  or  vg  the  phase 
difference  8g  can  be  calculated  by 

sin  8g  =  Corr  [ug,  H(vg)]  =  -Corr  [ vg ,  H(ug)],  (35) 

where  the  correlation  function  is  defined  as 


Corr  <“> 

Combining  (31)  and  (35)  gives  the  circular  polarization  param¬ 
eter  Q .  This  technique  of  using  the  covariance  of  u '  and  H(vr) 
for  calculating  Q  was  also  indicated  by  Eckermann  et  al.  [1994] 
for  the  total  perturbation  velocities  (u ' ,  v').  Such  an  approach 
leads  to  a  physically  meaningful  Q  only  when  there  exists  a 
dominant  wave  component,  as  indicated  by  Vincent  and  Fritts 
[1987]  by  a  condition  of  sufficiently  narrow  bands.  Here  our 
( ug ,  vg)  profiles  are  derived  from  the  EMD  method  that 
guarantees  narrow  bands.  Therefore  present  sin  8g  or  Q  de¬ 
rived  from  (35)  is  accurate.  The  phase  difference  8g  is  related 
to  the  gravity  wave  parameters  through  [Zhu,  1987] 


where  /  is  the  Coriolis  parameter,  co  is  the  Doppler-shifted 
frequency,  (k,  €)  is  the  horizontal  wavenumber  vector,  and 
K2  =  k2  +  i2.  A  positive  correlation  between  ug  and  H(vg) 
corresponds  to  upward  propagating  gravity  waves  in  the  north¬ 
ern  hemisphere.  Furthermore,  (37)  indicates  that  |sin  ^  1 
under  either  one  of  the  following  three  conditions:  /c  — >  0, 
i  —*•  0,  or  a)  — >  /.  For  a  fast  traveling  gravity  wave  with 
a)  —>  N  /,  then  |  sin  <sc  1  which  corresponds  to  a  linear 
polarized  wave  motion  due  to  the  weak  inertial  effect. 

Figure  11  shows  the  scattered  plots  of  Corr  [ug,  H(vg)]  and 
the  major  axis  orientation  rg,  and  we  see  that  the  upward 
propagating  gravity  waves  are  dominant  in  the  middle  atmo¬ 
sphere  which  is  consistent  with  many  previous  studies  [e.g., 
Vincent ,  1984;  Hirota  and  Niki,  1985].  Furthermore,  Figure  lib 
indicate  that  above  —50  km  there  is  an  anisotropy  in  the 
horizontal  directions  of  wave  propagation  with  an  NW/SE  bias. 
Such  an  NW/SE  bias  in  the  northern  hemisphere  winter  was 
also  detected  previously  by  Hamilton  [1991]. 

6.  Conclusions 

We  have  explored  the  recently  developed  Empirical  Mode 
Decomposition  (EMD)  method  and  applied  it  to  analyzing 
gravity  wave  characteristics  in  the  middle  atmosphere.  An  im¬ 
portant  connection  is  established  in  this  paper  between  the 
fundamental  Intrinsic  Mode  Functions  (IMFs)  derived  from 
the  EMD  method  and  WKB  solutions  of  a  dispersive- 
dissipative  wave  equation.  Since  the  WKB  solution  corre¬ 
sponds  to  a  wave  packet  whose  structure  and  wavenumber  are 
internally  determined  from  the  data,  the  IMF  representations 
of  wind  profiles  provide  useful  insights  into  physical  processes 
in  the  middle  atmosphere  where  the  dispersive-dissipative 
wave  phenomena  are  dominant.  The  positive  distribution  func¬ 
tion  P(zk ,  m€)  introduced  in  section  3  provides  a  quantitative 


16,560 


ZHU  ET  AL.:  GRAVITY  WAVE  CHARACTERISTICS  IN  THE  MIDDLE  ATMOSPHERE 


description  of  space-wavenumber  localization  of  a  data  series. 
Application  to  the  observed  wind  profiles  derived  from  the 
D  YANA  campaign  leads  to  an  identification  of  what  appear  to 
be  unphysical  small  scale  oscillations  by  falling  spheres  embed¬ 
ded  in  the  wind  profiles  above  60  km.  It  is  found  that  for  the 
vertical  wavenumber  m  <  3  km-1  ( Lz  >  2  km)  the  previously 
observed  m~3  Fourier  spectra  could  be  a  manifestation  of  a 
linear  wave  packet  whose  characteristic  vertical  wavenumber 
decreases  with  altitude.  For  small  vertical  scale  disturbances 
with  m  >  3  km-1  ( Lz  <  2  km)  a  near  -3  slope  in  the 
marginal  distribution  exists  locally  in  the  middle  atmosphere 
with  a  great  degree  of  universality,  suggesting  that  nonlinear 
energy  cascade  processes  might  dominate  the  spectral  forma¬ 
tion.  It  is  also  shown  that  the  hodographic  analysis  by  the 
Stokes  parameters  advocated  by  Vincent  and  Fritts  [1987]  and 
Eckermann  and  Vincent  [1989]  can  also  be  applied  to  the  IMFs 
derived  from  the  EMD  method. 
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